Tarea series de tiempo - Carlos Vesga, Santigo Buitrago

Contents

Tarea series de tiempo - Carlos Vesga, Santigo Buitrago#

1. Considere la serie de tiempo asociada con los futuros de la criptomoneda Bitcoin desde que comenzó a comercializarse hasta la fecha del día de hoy. Utilice la API de Yahoo Finance para obtener esta serie de tiempo.#

Lectura de datos y EDA#

import warnings
warnings.filterwarnings("ignore")
import yfinance as yf
from yahoo_fin.stock_info import get_data
from datetime import datetime, timedelta
import seaborn as sns
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from numpy import log
import numpy as np, pandas as pd
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
Warning - Certain functionality 
             requires requests_html, which is not installed.
             
             Install using: 
             pip install requests_html
             
             After installation, you may have to restart your Python session.
msft = yf.Ticker("BTC-USD")
msft.info
{'name': 'Bitcoin',
 'startDate': 1278979200,
 'description': 'Bitcoin (BTC) is a cryptocurrency launched in 2010. Users are able to generate BTC through the process of mining. Bitcoin has a current supply of 19,683,106. The last known price of Bitcoin is 66,797.9989346 USD and is up 3.99 over the last 24 hours. It is currently trading on 11005 active market(s) with $43,228,757,125.25 traded over the last 24 hours. More information can be found at https://bitcoin.org/.',
 'maxAge': 86400,
 'priceHint': 2,
 'previousClose': 63441.96,
 'open': 63441.96,
 'dayLow': 62850.684,
 'dayHigh': 63825.55,
 'regularMarketPreviousClose': 63441.96,
 'regularMarketOpen': 63441.96,
 'regularMarketDayLow': 62850.684,
 'regularMarketDayHigh': 63825.55,
 'volume': 42898325504,
 'regularMarketVolume': 42898325504,
 'averageVolume': 34966998847,
 'averageVolume10days': 36316462166,
 'averageDailyVolume10Day': 36316462166,
 'marketCap': 1251204005888,
 'fiftyTwoWeekLow': 24797.168,
 'fiftyTwoWeekHigh': 73750.07,
 'fiftyDayAverage': 66794.59,
 'twoHundredDayAverage': 46248.73,
 'currency': 'USD',
 'fromCurrency': 'BTC',
 'toCurrency': 'USD=X',
 'lastMarket': 'CoinMarketCap',
 'coinMarketCapLink': 'https://coinmarketcap.com/currencies/bitcoin',
 'volume24Hr': 42898325504,
 'volumeAllCurrencies': 42898325504,
 'circulatingSupply': 19683776,
 'exchange': 'CCC',
 'quoteType': 'CRYPTOCURRENCY',
 'symbol': 'BTC-USD',
 'underlyingSymbol': 'BTC-USD',
 'shortName': 'Bitcoin USD',
 'longName': 'Bitcoin USD',
 'firstTradeDateEpochUtc': 1410912000,
 'timeZoneFullName': 'UTC',
 'timeZoneShortName': 'UTC',
 'uuid': '74397779-1589-3270-8c45-b7f1a7345b3a',
 'messageBoardId': 'finmb_BTC_CCC',
 'trailingPegRatio': None}

Podemos ver los datos basicos de la accion, tales como su valor actual, la fecha en que inicio a ser comercializada entre otros.

sns.set_theme()
sns.set_context("paper")
stock = 'BTC-USD'
resolution = '1d'
end_date = datetime.now()
start_date = 1278979200

Como queremos analizar la informacion desde el inicio de la serie de tiempo hasta el momento actual el end_date debe ser igual a “hoy” y el start_date debe coincidir con le fecha de inicio mostrada anteriormente. Como estamos interesados en predecir el valor del día nuestra resolucion es de un día.

def date_format(date_h):
    return date_h.strftime('%d/%m/%Y')
BTC_df = get_data(stock, start_date=start_date, end_date=end_date, interval=resolution, index_as_date=False)
BTC_df.head()
date open high low close adjclose volume ticker
0 2014-09-17 465.864014 468.174011 452.421997 457.334015 457.334015 21056800.0 BTC-USD
1 2014-09-18 456.859985 456.859985 413.104004 424.440002 424.440002 34483200.0 BTC-USD
2 2014-09-19 424.102997 427.834991 384.532013 394.795990 394.795990 37919700.0 BTC-USD
3 2014-09-20 394.673004 423.295990 389.882996 408.903992 408.903992 36863600.0 BTC-USD
4 2014-09-21 408.084991 412.425995 393.181000 398.821014 398.821014 26580100.0 BTC-USD
sns.lineplot(data=BTC_df, x=BTC_df.date, y=BTC_df.close);
_images/3c0f846ccaf3568d348569673ae911f42714ab614f0dafbc4c70976d3394242e.png

Al leer y graficar los datos vemos una gran volatilidad en el valor de los BTC, esto es esperado y cabe señalar que entorpece las predicciones realizadas.

import plotly.graph_objects as go

fig = go.Figure(data=[go.Candlestick(x = BTC_df.date,
                                     open = BTC_df.open, 
                                     high = BTC_df.high,
                                     low = BTC_df.low, 
                                     close = BTC_df.close)
                     ])
fig.update_layout(
    title="Bitcoin. (BTC)",
    xaxis_title="Day",
    yaxis_title="BTC-USD",
    font=dict(
        family="Courier New, monospace",
        size=12,
        color="RebeccaPurple"
    )
)
fig.update_layout(xaxis_rangeslider_visible=False)
fig

Se grafica tambien un candlestick para facilitar la visualizacion.

Pruebas estadisticas#

Antes de comenzar a pronosticar decidimos realizar algunas pruebas estadisticas para confirmar la viabilidad del modelo ARIMA.

result = adfuller(BTC_df.close)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: nan
p-value: nan

Al p-value ser mayor a 0.05 no rechazamos la hipotesis nula, por lo cual la serie estadisticamente no es estacionaria.

import plotly.graph_objs as go
from plotly.subplots import make_subplots
import numpy as np
from statsmodels.tsa.stattools import acf

def plotly_acf_with_confidence(data, lags):
    acf_vals, conf_int = acf(data, nlags=lags, alpha=0.05, fft=True)  # Obtener ACF y intervalos de confianza
    lower_bound = conf_int[:, 0]  # límite inferior
    upper_bound = conf_int[:, 1]  # límite superior
    
    traces = [
        go.Scatter(
            x=list(range(lags+1)), y=acf_vals, mode='lines', name='ACF',
            line=dict(color='blue'),
            fill='tozeroy',  # Relleno hasta el eje x
            fillcolor='rgba(0, 100, 255, 0.2)'
        ),
        go.Scatter(
            x=list(range(lags+1)) + list(range(lags+1))[::-1],  # x, luego x reversa
            y=np.concatenate([upper_bound, lower_bound[::-1]]),  # upper, luego lower reversa
            fill='toself', fillcolor='rgba(0, 100, 255, 0.3)', line=dict(color='rgba(255,255,255,0)'),
            showlegend=False  # Ocultar leyenda para este trazo
        )
    ]
    return traces

# Suponiendo que tienes 'BTC_df' cargado y listo, y la columna 'close' tiene los precios de cierre
fig = make_subplots(rows=3, cols=2, subplot_titles=("Original Series", "ACF of Original",
                                                    "1st Order Differencing", "ACF of 1st Differencing",
                                                    "2nd Order Differencing", "ACF of 2nd Differencing"))

# Añadiendo las series y los ACF con intervalos de confianza
fig.add_trace(go.Scatter(x=BTC_df.index, y=BTC_df.close, mode='lines', showlegend=False), row=1, col=1)
fig.add_traces(plotly_acf_with_confidence(BTC_df.close, 3497), rows=[1,1], cols=[2,2])

# Primera diferenciación
first_diff = BTC_df.close.diff().dropna()
fig.add_trace(go.Scatter(x=BTC_df.index[1:], y=first_diff, mode='lines', showlegend=False), row=2, col=1)
fig.add_traces(plotly_acf_with_confidence(first_diff, 3496), rows=[2,2], cols=[2,2])

# Segunda diferenciación
second_diff = BTC_df.close.diff().diff().dropna()
fig.add_trace(go.Scatter(x=BTC_df.index[2:], y=second_diff, mode='lines', showlegend=False), row=3, col=1)
fig.add_traces(plotly_acf_with_confidence(second_diff, 3495), rows=[3,3], cols=[2,2])

# Actualizar layout para mejorar la visualización
fig.update_layout(height=1200, width=1000, title_text="BTC Price Analysis with Interactive ACF Plots")
fig.show()

Graficamos las auto-correlaciones de la grafica original, la de difereciacion de orden 1 y de orden 2. El decaimiento geometrico en la primera grafica nos hace confirmar que la serie de tiempo no es estacionaria, ademas podemos ver que la grafica de diferenciacion de orden 2 entra al lado negativo con rapidez, esto puede ser debido a que esta sobrediferenciada.

2. Repita TODOS los pasos indicados en esta sección para encontrar modelos ARIMA para predecir el precio de Bitcoin con los siguientes horizontes: 7, 14, 21 y 28 días. Utilizar siempre predicciones usando rolling con ventana de predicción continua de un día. Cualquier cantidad de pasos extra para enriquecer su análisis predictivo serán aceptados siempre y cuando sean acordes con lo que indica la teoría de análisis de series de tiempo.#

from statsmodels.tsa.arima.model import ARIMA
import warnings
from statsmodels.graphics.tsaplots import plot_predict
warnings.filterwarnings("ignore")

Hacemos el llamado de la libreria.

def arima_rolling(history, test, orderb):
    
    predictions = list()
    for t in range(len(test)):
        model = ARIMA(history, order=orderb)
        model_fit = model.fit()
        output = model_fit.forecast()
        yhat = output[0]
        predictions.append(yhat)
        obs = test[t]
        history.append(obs)
        
    return predictions

Definimos la funcion para trabajar en la metodologia rolling forecast.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score

def forecast_accuracy(predictions_list, actual, title_prefix=""):
    # Asumiendo que 'actual' es una serie o arreglo que contiene los valores reales
    
    results = []
    fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(10, 18))
    metrics = ['MAE', 'MSE', 'MAPE', 'RMSE', 'R2']
    colors = ['blue', 'green', 'red', 'purple', 'orange']
    
    for idx, forecast in enumerate(predictions_list):
        days_predicted = len(forecast)
        str_name = f"{title_prefix} {days_predicted} days"
        
        # Calculo de errores
        mape = np.mean(np.abs(forecast - actual[:days_predicted]) / np.abs(actual[:days_predicted]))  # MAPE
        mae = np.mean(np.abs(forecast - actual[:days_predicted]))  # MAE
        rmse = np.mean((forecast - actual[:days_predicted])**2)**.5  # RMSE
        mse = np.mean((forecast - actual[:days_predicted])**2)  # MSE
        r2 = r2_score(actual[:days_predicted], forecast)  # R2
        
        results.append({
            'Days': days_predicted,
            'MAE': mae,
            'MSE': mse,
            'MAPE': mape,
            'RMSE': rmse,
            'R2': r2
        })
    
    # Crear un DataFrame para los resultados
    df_results = pd.DataFrame(results)
    
    # Graficar cada métrica
    for i, metric in enumerate(metrics):
        axes[i].plot(df_results['Days'], df_results[metric], label=metric, color=colors[i], marker='o')
        axes[i].set_title(f"{metric} by days predicted")
        axes[i].set_xlabel('Days')
        axes[i].set_ylabel(metric)
        axes[i].legend()
    
    plt.tight_layout()
    plt.show()
    
    return df_results
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[14], line 4
      2 import pandas as pd
      3 import matplotlib.pyplot as plt
----> 4 from sklearn.metrics import r2_score
      6 def forecast_accuracy(predictions_list, actual, title_prefix=""):
      7     # Asumiendo que 'actual' es una serie o arreglo que contiene los valores reales
      9     results = []

ModuleNotFoundError: No module named 'sklearn'

Se define esta funcion para el calculo de las metricas de error. Modificamos la funcion enseñada en clase para que acepte listas y grafique cada metrica.

import numpy as np
import scipy.stats as stats
from statsmodels.stats.diagnostic import acorr_ljungbox

def analyze_residuals(predictions_list, actual_data):
    results = []

    for idx, prediction in enumerate(predictions_list):
        # Asumimos que 'actual_data' tiene suficientes datos para comparar con cada 'prediction'
        actual_subset = actual_data[:len(prediction)]
        
        # Calcular residuos
        residuals = np.array(actual_subset) - np.array(prediction)
        
        # Prueba de normalidad (Shapiro-Wilk)
        shapiro_test = stats.shapiro(residuals)
        shapiro_pvalue = shapiro_test.pvalue
        
        # Prueba de independencia (Ljung-Box)
        # El número de lags recomendado por defecto es min(10, len(residuals) // 5)
        lb_lags = min(10, len(residuals) // 5)
        lb_test = acorr_ljungbox(residuals, lags=[lb_lags], return_df=True)
        
        # Almacenar resultados
        results.append({
            'Prediction Length': len(prediction),
            'Shapiro P-Value': shapiro_pvalue,
            'Ljung-Box P-Value': lb_test['lb_pvalue'].iloc[0]
        })
        
        # Imprimir conclusiones
        print(f"Residuos para la predicción de {len(prediction)} días:")
        print(f"Test de Shapiro-Wilk P-Value: {shapiro_pvalue:.4f}")
        if shapiro_pvalue < 0.05:
            print("Rechazar la hipótesis nula: Los residuos no son normales.")
        else:
            print("No rechazar la hipótesis nula: Los residuos parecen ser normales.")
        
        print(f"Test de Ljung-Box P-Value: {lb_test['lb_pvalue'].iloc[0]:.4f}")
        if lb_test['lb_pvalue'].iloc[0] < 0.05:
            print("Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.")
        else:
            print("No rechazar la hipótesis nula: Los residuos parecen ser independientes.")
        print("\n")
        
    return results

Se crea esta funcion para facilitar las pruebas de normalidad e independencia de los residuales de cada intervalo de prediccion.

n_BTC = len(BTC_df.close); n_test = 28 
train_size = n_BTC - n_test
train = BTC_df.close[:train_size]
dates_train = BTC_df.date[:train_size]
test_4w = BTC_df.close[train_size:train_size + n_test] 
dates_4w = BTC_df.date[train_size:train_size + n_test] 
print("train:", train.shape)
print("test_4w:", test_4w.shape)
train: (3470,)
test_4w: (28,)

Definimos y dividimos el training y el test

train_df = BTC_df[["close"]][:train_size]
test_4w_df = BTC_df[["close"]][train_size:train_size + n_test] 

Organizamos el df en train y test dejando los ultimos 28 días para el testing.

best_aic = np.inf

best_order = None
best_mdl = None

pq_rng = range(5)
d_rng  = range(3)

for i in pq_rng:
    for d in d_rng:
        for j in pq_rng:
            try:
                # print(i, d, j)
                tmp_mdl = ARIMA(train, order=(i,d,j)).fit()
                tmp_aic = tmp_mdl.aic
                if tmp_aic < best_aic:
                    best_aic = tmp_aic
                    best_order = (i, d, j)
                    best_mdl = tmp_mdl
            except: continue

Hacemos la hyperparametrizacion del modelo con los parametros moviendose del 1 al 10 usando la metrica aic.

print('aic: {:6.5f} | order: {}'.format(best_aic, best_order))
aic: 56295.88868 | order: (3, 2, 4)

Podemos ver que los mejores parametros segun la metrica aic son (3, 2, 4). Esto es un poco contradictorio con el analisis hecho con las graficas de autocorrelaciones pero trabajaremos con este orden y segun las pruebas estadisticas decidiremos si bajar la diferenciacion a una o mantenerla.

model_aic = ARIMA(train, order=best_order)
model_fit_aic = model_aic.fit()

Entrenamos el modelo con los mejores parametros encontrados.

plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit_aic, 1, ax=ax);
plt.show();
_images/91e25f78a275f5d7f64e7e5968c933fcf032dc769002050f7fbb5afab4d04129.png

Graficamos el intervalo de confianza del 95% para nuestras predicciones.

test_range = range(4)
predictions_roll_aic = []
for i in test_range:
    test = test_4w.head(7*(i+1))
    test = test.tolist()
    yhat  = arima_rolling(train.tolist(), test, best_order)
    predictions_roll_aic.append(yhat)

En esta funcion estamos haciendo el rolling forecast para los intervalos de 7,14,21 y 28 días, nos ayudamos en el hecho de ser todos multiplos del 7. Guardamos los resultados de cada pronostico en una lista.

results = analyze_residuals(predictions_roll_aic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.7373
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1296
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.4232
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.3658
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.3687
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5013
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.7656
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5265
No rechazar la hipótesis nula: Los residuos parecen ser independientes.

Al utilizar la funcion definida anteriormente podemos ver que todos los supuestos se cumplen para cada residual, por lo que concluimos que el modelo capta la variabilidad de los datos y es posible confiar en las predicciones.

3. Repita el paso 2 ahora sin utilizar rolling. Esto es, realice el pronóstico solo utilizando forecast() para los diferentes horizontes de predicción, 7, 14, 21 y 28 días.#

test_range = range(4)
predictions_for_aic = []

for i in test_range:
    # Hacer las predicciones para el siguiente intervalo de tiempo
    num_steps = 7 * (i + 1)  # Calculando el número de pasos a predecir
    yhat = model_fit_aic.forecast(steps=num_steps)
    predictions_for_aic.append(yhat)
results = analyze_residuals(predictions_for_aic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.3997
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0792
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.1172
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0402
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.


Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.1362
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1099
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.0445
Rechazar la hipótesis nula: Los residuos no son normales.
Test de Ljung-Box P-Value: 0.0416
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.

Observamos que a diferencia del rolling forecast el supuesto de independencia no se cumple para algunos intervalos de tiempo, se necesita una investigacion en el tema pero para ambitos educativos trabajaremos asi.

4. Realice tablas de error para los ítems 1 y 2, utilizando las métricas: MAPE, MAE, RMSE, MSE, R2. Además, agregue el gráfico de correlación entre la observación real y su predicción en el test.#

Rolling Forecast#

roll_aic = forecast_accuracy(predictions_roll_aic, test_4w, title_prefix="Forecast")
_images/90c2f05f0df623fbb64c6304e508784769b98a01589825ec85101f6c48b7c433.png

Graficamos y guardamos los resultados conseguidos, por extraño que pareza parece que los resultados mejoran al aumentar el intervalo de prediccion.

Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:

  • Datos limitados

  • Sobreajuste

  • Variacion Aleatoria

import plotly.graph_objs as go
from plotly.subplots import make_subplots
import numpy as np

def plot_prediction_correlation(actual, predictions_list):
    # Crear subplots, uno para cada conjunto de predicciones
    num_predictions = len(predictions_list)
    fig = make_subplots(rows=num_predictions, cols=1, 
                        subplot_titles=[f'Correlacion de {len(pred)} dias predichos' for pred in predictions_list],
                        vertical_spacing=0.1)

    for i, predictions in enumerate(predictions_list):
        # Asegurarse de que los datos reales y las predicciones estén alineados
        actual_subset = actual[:len(predictions)]

        # Calcular coeficientes de correlación
        correlation_matrix = np.corrcoef(actual_subset, predictions)
        correlation_xy = correlation_matrix[0,1]
        r_squared = correlation_xy**2

        # Añadir gráfico de dispersión para cada conjunto de predicciones
        fig.add_trace(
            go.Scatter(
                x=actual_subset, y=predictions, mode='markers',
                marker=dict(size=8, color='blue', opacity=0.5),
                name=f'R^2: {r_squared:.2f}'
            ), 
            row=i+1, col=1
        )
        
        # Ajustar una línea de regresión lineal
        z = np.polyfit(actual_subset, predictions, 1)
        p = np.poly1d(z)
        fig.add_trace(
            go.Scatter(
                x=actual_subset, y=p(actual_subset), mode='lines',
                line=dict(color='red', dash='dash'),
                name='Fit'
            ),
            row=i+1, col=1
        )

        # Actualizar ejes para cada subplot
        fig.update_xaxes(title_text="Valores Reales", row=i+1, col=1)
        fig.update_yaxes(title_text="Valores Predichos", row=i+1, col=1)

    # Actualizar layout del gráfico
    fig.update_layout(
        height=300 * num_predictions, width=800,
        title_text="Correlacion entre predichos y verdaderos",
        showlegend=False
    )
    
    fig.show()

Creamos una funcion para graficar las correlaciones entre el valor real y predicho de cada intervalo de tiempo

plot_prediction_correlation(test_4w, predictions_roll_aic)

Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, los mejores R^2 estan para el intervalo de 28 y 14 dias.

Forecast#

for_aic = forecast_accuracy(predictions_for_aic, test_4w, title_prefix="Forecast")
_images/aa95f512e7fc6101b915f156052940c4acbb9f0a2d0dee8302fbac99156e59fb.png

Parece que el mejor intervalo es a 21 días.

Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:

  • Datos limitados

  • Sobreajuste

  • Variacion Aleatoria

plot_prediction_correlation(test_4w, predictions_for_aic)

Todos los intervalos tienen una tendencia positiva,y el mejor R^2 esta a 14 días.

Comparaciones#

import matplotlib.pyplot as plt

def plot_comparative_metrics(df1, df2):
    metrics = ['MAE', 'MSE', 'MAPE', 'RMSE', 'R2']
    num_metrics = len(metrics)
    intervals = df1.index  # Asumimos que el índice es el intervalo de días
    width = 0.35  # Ancho de las barras

    # Crear una figura grande para contener todos los subplots
    fig, axes = plt.subplots(nrows=num_metrics, ncols=1, figsize=(10, 6 * num_metrics))

    for i, metric in enumerate(metrics):
        ax = axes[i]  # Acceder al subplot correspondiente
        # Graficar las barras para DF1
        ax.bar(intervals - width/2, df1[metric], width, label='Rolling', color='blue')
        # Graficar las barras para DF2
        ax.bar(intervals + width/2, df2[metric], width, label='Forecast', color='orange')

        ax.set_title(f'Comparación de {metric}')
        ax.set_ylabel(metric)
        ax.legend()

    axes[-1].set_xlabel('Intervalos de días')  # Etiqueta del eje x solo en el último gráfico
    axes[-1].set_xticks(intervals)  # Asegurar que todos los intervalos estén etiquetados

    plt.tight_layout()
    plt.show()
plot_comparative_metrics(roll_aic, for_aic)
_images/1b94a9b583eabc3d1a17ee40e0aa8b1385b250644815a25a1ba6480f7aa94723.png

En terminos generales el rolling forecast es mejor, debido a que minimiza la mayoria de los errores.

5. Repita el análisis desarrollado en los pasos anteriores, considerando ahora el criterio de inferencia Bayesiana (BIC) y el criterio de información de Hannan–Quinn (HQIC) para encontrar el mejor modelo ARIMA y, compare los errores con aquellos obtenidos con el criterio de Akaike.#

BIC#

best_bic = np.inf

best_order = None
best_mdl = None

pq_rng = range(5)
d_rng  = range(3)

for i in pq_rng:
    for d in d_rng:
        for j in pq_rng:
            try:
                # print(i, d, j)
                tmp_mdl = ARIMA(train, order=(i,d,j)).fit()
                tmp_bic = tmp_mdl.bic
                if tmp_bic < best_bic:
                    best_bic = tmp_bic
                    best_order = (i, d, j)
                    best_mdl = tmp_mdl
            except: continue

Volvemos a correr la hyperparametrizacion esta vez optimizando el bic

print('bic: {:6.5f} | order: {}'.format(best_bic, best_order))
bic: 56327.16155 | order: (0, 2, 1)
model_bic = ARIMA(train, order=best_order)
model_fit_bic = model_bic.fit()
plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit_bic, 1, ax=ax);
plt.show();
_images/60c133405e4fd5528902e9014e89a838a4cf0c2e5018266c1dac3b639f3220b6.png

Graficamos el intervalo de confianza del 95% para nuestro nuevo modelo.

Rolling Forecast#

test_range = range(4)
predictions_roll_bic = []
for i in test_range:
    test = test_4w.head(7*(i+1))
    test = test.tolist()
    yhat  = arima_rolling(train.tolist(), test, best_order)
    predictions_roll_bic.append(yhat)

Hacemos el rolling forecast para el nuevo modelo

results = analyze_residuals(predictions_roll_bic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.9388
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1408
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.8911
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.3486
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.8613
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5638
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.9871
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.6047
No rechazar la hipótesis nula: Los residuos parecen ser independientes.

Cumple todos los supuestos, por lo que concluimos que el modelo capta la variabilidad de los datos y es posible confiar en las predicciones.

roll_bic = forecast_accuracy(predictions_roll_bic, test_4w, title_prefix="Forecast")
_images/3dae36f64ac7566b2eb36a1e4aa48307b82d4bd1a4e654a53b3f028cb48946d7.png

Graficamos y guardamos los resultados conseguidos, por extraño que pareza parece que los resultados mejoran al aumentar el intervalo de prediccion.

Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:

  • Datos limitados

  • Sobreajuste

  • Variacion Aleatoria

plot_prediction_correlation(test_4w, predictions_roll_bic)

Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 28 dias.

Forecast#

test_range = range(4)
predictions_for_bic = []

for i in test_range:
    # Hacer las predicciones para el siguiente intervalo de tiempo
    num_steps = 7 * (i + 1)  # Calculando el número de pasos a predecir
    yhat = model_fit_bic.forecast(steps=num_steps)
    predictions_for_bic.append(yhat)
results = analyze_residuals(predictions_for_bic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.4822
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0864
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.1098
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0435
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.


Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.1307
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1106
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.0391
Rechazar la hipótesis nula: Los residuos no son normales.
Test de Ljung-Box P-Value: 0.0419
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.

Observamos que a diferencia del rolling forecast el supuesto de independencia no se cumple para algunos intervalos de tiempo, se necesita una investigacion en el tema pero para ambitos educativos trabajaremos asi.

for_bic = forecast_accuracy(predictions_for_bic, test_4w, title_prefix="Forecast")
_images/c49edb9ed769b089dca30e6cb0c6351f2709c2426c137583c4796971310e0467.png

Parece que el mejor intervalo es a 21 días.

Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:

  • Datos limitados

  • Sobreajuste

  • Variacion Aleatoria

plot_prediction_correlation(test_4w, predictions_for_bic)

Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 14 dias.

Comparaciones#

plot_comparative_metrics(roll_bic, for_bic)
_images/974d3dfeab5e97d4156bcfe937db0356e37f6fbb4c4688d050ccdafe7227c972.png

En terminos generales el rolling forecast es mejor, debido a que minimiza la mayoria de los errores.

HQIC#

best_hqic = np.inf

best_order = None
best_mdl = None

pq_rng = range(5)
d_rng  = range(3)

for i in pq_rng:
    for d in d_rng:
        for j in pq_rng:
            try:
                # print(i, d, j)
                tmp_mdl = ARIMA(train, order=(i,d,j)).fit()
                tmp_hqic = tmp_mdl.hqic
                if tmp_hqic < best_hqic:
                    best_hqic = tmp_hqic
                    best_order = (i, d, j)
                    best_mdl = tmp_mdl
            except: continue

Hyperparametrizamos optimizando el HQIC

print('HQIC: {:6.5f} | order: {}'.format(best_hqic, best_order))
HQIC: 56313.45958 | order: (3, 2, 4)
model_hqic = ARIMA(train, order=best_order)
model_fit_hqic = model_hqic.fit()
plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit_hqic, 1, ax=ax);
plt.show();
_images/91e25f78a275f5d7f64e7e5968c933fcf032dc769002050f7fbb5afab4d04129.png

Graficamos el intervalo de confianza del 95%

Rolling Forecast#

test_range = range(4)
predictions_roll_hqic = []
for i in test_range:
    test = test_4w.head(7*(i+1))
    test = test.tolist()
    yhat  = arima_rolling(train.tolist(), test, best_order)
    predictions_roll_hqic.append(yhat)

Hacemos el rolling forecast para el nuevo modelo

results = analyze_residuals(predictions_roll_hqic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.7373
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1296
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.4232
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.3658
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.3687
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5013
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.7656
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.5265
No rechazar la hipótesis nula: Los residuos parecen ser independientes.

Cumple todos los supuestos, por lo que concluimos que el modelo capta la variabilidad de los datos y es posible confiar en las predicciones.

roll_hqic = forecast_accuracy(predictions_roll_hqic, test_4w, title_prefix="Forecast")
_images/90c2f05f0df623fbb64c6304e508784769b98a01589825ec85101f6c48b7c433.png

Graficamos y guardamos los resultados conseguidos, por extraño que pareza parece que los resultados mejoran al aumentar el intervalo de prediccion.

Vemos que el R^2 tiene valores negativos, esto ocurre debido a varios factores como:

  • Datos limitados

  • Sobreajuste

  • Variacion Aleatoria

plot_prediction_correlation(test_4w, predictions_roll_hqic)

Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 28 dias.

Forecast#

test_range = range(4)
predictions_for_hqic = []

for i in test_range:
    # Hacer las predicciones para el siguiente intervalo de tiempo
    num_steps = 7 * (i + 1)  # Calculando el número de pasos a predecir
    yhat = model_fit_hqic.forecast(steps=num_steps)
    predictions_for_hqic.append(yhat)
results = analyze_residuals(predictions_for_hqic,test_4w)
Residuos para la predicción de 7 días:
Test de Shapiro-Wilk P-Value: 0.3997
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0792
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 14 días:
Test de Shapiro-Wilk P-Value: 0.1172
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.0402
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.


Residuos para la predicción de 21 días:
Test de Shapiro-Wilk P-Value: 0.1362
No rechazar la hipótesis nula: Los residuos parecen ser normales.
Test de Ljung-Box P-Value: 0.1099
No rechazar la hipótesis nula: Los residuos parecen ser independientes.


Residuos para la predicción de 28 días:
Test de Shapiro-Wilk P-Value: 0.0445
Rechazar la hipótesis nula: Los residuos no son normales.
Test de Ljung-Box P-Value: 0.0416
Rechazar la hipótesis nula: Los residuos muestran autocorrelación significativa.

Observamos que a diferencia del rolling forecast el supuesto de independencia no se cumple para algunos intervalos de tiempo, se necesita una investigacion en el tema pero para ambitos educativos trabajaremos asi.

for_hqic = forecast_accuracy(predictions_for_hqic, test_4w, title_prefix="Forecast")
_images/aa95f512e7fc6101b915f156052940c4acbb9f0a2d0dee8302fbac99156e59fb.png

Parece que el mejor intervalo es a 21 días.

plot_prediction_correlation(test_4w, predictions_for_hqic)

Vemos que a excepcion de el intervalo de 7 días todoso tienen una tendencia positiva, ademas, el mejor R^2 esta en el intervalo de 14 dias.

Comparaciones#

plot_comparative_metrics(roll_aic, for_aic)
_images/1b94a9b583eabc3d1a17ee40e0aa8b1385b250644815a25a1ba6480f7aa94723.png

En terminos generales el rolling forecast es mejor, debido a que minimiza la mayoria de los errores.

Comparaciones#

import matplotlib.pyplot as plt

def plot_comparative_metrics_3(df1, df2, df3):
    metrics = ['MAE', 'MSE', 'MAPE', 'RMSE', 'R2']
    num_metrics = len(metrics)
    intervals = df1.index  # Asumimos que el índice es el intervalo de días
    width = 0.20  # Ancho de las barras, ajustado para acomodar tres modelos

    # Crear una figura grande para contener todos los subplots
    fig, axes = plt.subplots(nrows=num_metrics, ncols=1, figsize=(10, 6 * num_metrics))

    for i, metric in enumerate(metrics):
        ax = axes[i]  # Acceder al subplot correspondiente
        # Graficar las barras para DF1
        ax.bar(intervals - width, df1[metric], width, label='AIC', color='blue')
        # Graficar las barras para DF2
        ax.bar(intervals, df2[metric], width, label='BIC', color='orange')
        # Graficar las barras para DF3
        ax.bar(intervals + width, df3[metric], width, label='HQIC', color='green')

        ax.set_title(f'Comparación de {metric}')
        ax.set_ylabel(metric)
        ax.legend()

    axes[-1].set_xlabel('Intervalos de días')  # Etiqueta del eje x solo en el último gráfico
    axes[-1].set_xticks(intervals)  # Asegurar que todos los intervalos estén etiquetados

    plt.tight_layout()
    plt.show()
plot_comparative_metrics_3(roll_aic, roll_bic, roll_hqic)
_images/ec12a91fb1e8e1fc2fb412de8daec28dfca5f56ca035670fa168c97bebba1a53.png